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Abstract 

In vitro cultures of endothelial cells are a widely used model system of the collective behavior of endothelial cells during 
vasculogenesis and angiogenesis. When seeded in an extracellular matrix, endothelial cells can form blood vessel-like 
structures, including vascular networks and sprouts. Endothelial morphogenesis depends on a large number of chemical 
and mechanical factors, including the compliancy of the extracellular matrix, the available growth factors, the adhesion of 
cells to the extracellular matrix, cell-cell signaling, etc. Although various computational models have been proposed to 
explain the role of each of these biochemical and biomechanical effects, the understanding of the mechanisms underlying 
in vitro angiogenesis is still incomplete. Most explanations focus on predicting the whole vascular network or sprout from 
the underlying cell behavior, and do not check if the same model also correctly captures the intermediate scale: the pairwise 
cell-cell interactions or single cell responses to ECM mechanics. Here we show, using a hybrid cellular Potts and finite 
element computational model, that a single set of biologically plausible rules describing (a) the contractile forces that 
endothelial cells exert on the ECM, (b) the resulting strains in the extracellular matrix, and (c) the cellular response to the 
strains, suffices for reproducing the behavior of individual endothelial cells and the interactions of endothelial cell pairs in 
compliant matrices. With the same set of rules, the model also reproduces network formation from scattered cells, and 
sprouting from endothelial spheroids. Combining the present mechanical model with aspects of previously proposed 
mechanical and chemical models may lead to a more complete understanding of in vitro angiogenesis. 
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Introduction 

How the behavior of cells in a multicellular organism is 
coordinated to form structured tissues, organs and whole 
organisms, is a central question in developmental biology. Keys 
to answering this question are chemical and mechanical cell-cell 
communication and the biophysics of self- organization. Cells 
exchange information by means of diffusing molecular signals, and 
by membrane-bound molecular signals for which direct cell-cell 
contact is required. In general, these developmental signals are 
short-lived and move over short distances. The extracellular 
matrix (ECM), the jelly or hard materials that cells secrete, 
provides the micro-environment the cells live in. Apart from its 
supportive function, the ECM mediates molecular [1] and 
biomechanical [2] signals between cells. Mechanical signals, in 
the form of tissue strains and stresses to which cells respond [3] , 
can act over long distances and integrate mechanical information 
over the whole tissue [4], and also mediate short-range, 



mechanical cell-cell communication [2]. How such mechanical 
cell-cell communication via the ECM can coordinate the self- 
organization of cells into tissues is still poorly understood. Here we 
propose a cell-based model of endothelial cell motility on 
compliant matrices to address this problem. 

A widely used approach to study the role of cell-ECM 
interactions in coordinating collective cell behavior is to isolate 
cells (e.g., endothelial cells isolate from bovine aortae or from 
human umbilical cords or foreskins) and culture them on top of or 
inside an artificial or natural ECM (e.g., Matrigel). This makes it 
possible to study the intrinsic ability of cells to form tissues in 
absence of potential organizing signals or pre-patterns from 
adjacent tissues. A problem particularly well-studied in cell 
cultures is the ability of endothelial cells to form blood vessel-like 
structures, including the formation of vascular-like networks from 
dispersed cells and the sprouting of spheroids. To this end, cell 
cultures can be initialized with a dispersion of endothelial cells on 
top of an ECM material (e.g., Matrigel, collagen, or fibrin) [5,6], 
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Author Summary 

During the embryonic development of multicellular 
organisms, millions of cells cooperatively build structured 
tissues, organs and whole organisms, a process called 
morphogenesis. How the behavior of so many cells is 
coordinated to produce complex structures is still incom- 
pletely understood. Most biomedical research focuses on 
the molecular signals that cells exchange with one 
another. It has now become clear that cells also commu- 
nicate biomechanically during morphogenesis. In cell 
cultures, endothelial cells — the building blocks of blood 
vessels — can organize into structures resembling networks 
of capillaries. Experimental work has shown that the 
endothelial cells pull onto the protein gel that they live in, 
called the extracellular matrix. On sufficiently compliant 
matrices, the strains resulting from these cellular pulling 
forces slow down and reorient adjacent cells. Here we 
propose a new computational model to show that this 
simple form of mechanical cell-cell communication suffices 
for reproducing the formation of blood vessel-like struc- 
tures in cell cultures. These findings advance our under- 
standing of biomechanical signaling during morphogene- 
sis, and introduce a new set of computational tools for 
modeling mechanical interactions between cells and the 
extracellular matrix. 

with endothelial spheroids embedded within the ECM [7,8], or 
with confluent endothelial monolayers [9—11]. Although the 
conditions required for vascular-like development in these in vitro 
culture systems are well established, the mechanisms driving 
pattern formation of endothelial cells are heavily debated, and a 
wide range of plausible mechanisms has been proposed in the form 
of mathematical and computational models reproducing aspects of 
angiogenesis (reviewed in [12-14]). 

Typical ingredients of network formation models are (a) an 
attractive force between endothelial cells, which is (b) proportional 
to the cell density, and (c) inhibited or attenuated at higher cellular 
densities. The attractive force can be due to mechanical traction or 
due to chemotaxis. Manoussaki, Murray, and coworkers [15,16] 
proposed a mechanical model of angiogenic network formation, 
based on the Oster and Murray [17,18] continuum mechanics 
theory of morphogenesis. In their model, endothelial cells exert a 
uniform traction force on the ECM, dragging the ECM and the 
associated endothelial cells towards them. The traction forces 
saturated at a maximum cell density. Namy and coworkers [19] 
replaced the endothelial cells' passive motion along with the ECM 
for active cell motility via haptotaxis, in which cells move actively 
towards higher concentrations of the ECM. Both models also 
included a strain-biased random walk term for the endothelial 
cells, but they found that it had little effect on network formation; 
the mechanism was dominated by cell aggregation. In their model 
based on chemotaxis, Preziosi and coworkers [20,21] assumed that 
cells attract one another via the secreted chemoattractant VEGF. 
Due to diffusion and first-order degradation, the chemoattractant 
forms exponential gradients around cells leading to cell aggrega- 
tion in much the same way as that assumed in the Manoussaki and 
Namy models. These chemotaxis-based hypotheses formed the 
basis for a series of cell-based models based on the cellular Potts 
model (CPM). Assuming chemotactic cell-cell attraction, and a 
biologically-plausible overdamped cell motility, the cells in these 
CPM models form round aggregates, in accordance with the 
Keller-Segel model of cell aggregation [22]. Additional assump- 
tions, including an elongated cell shape [23] or contact inhibition 
of chemotaxis [24] are needed to transform these circular 



aggregates into vascular-like network patterns. Related network 
formation models studied the role of ECM-bound growth factors 
[25—27] and a range of additional secreted and exogenous growth 
factors [27], and studied the ability of the contact-inhibition 
mechanism to produce three-dimensional blood-vessel-like struc- 
tures [28]. Szabo and coworkers found that in culture, astroglia- 
related rat C6 cells and muscle-related mouse C2C12 cells 
organize into network-like structures on rigid culture substrates 
[29], such that ECM-density or chemoattractant gradients are 
excluded. They proposed a model where cells were preferentially 
attracted to or preferentially adhered to locally elongated 
structures. As an alternative mechanism for "gel-free" network 
formation it was found that elongated cells can also produce 
networks in absence of chemoattractant gradients [30]. 

Paradoxically, despite the diverse assumptions underlying the 
mathematical models proposed for vascular network formation, 
many are at least partly supported by experimental evidence. This 
suggests that a combination of chemotaxis, and chemical and 
mechanical cell-ECM interactions drives network formation, or 
that each alternative mechanism operates in a different tissue, 
developmental stage, or culture condition. A problem is that one 
mathematical representation may represent a range of equivalent 
alternative underlying mechanisms. For example, a model 
representing cell-cell attraction cannot distinguish between che- 
motaxis-based cellular attraction [20,21,23,24], attraction via 
haptotaxis [19], direct mechanical attraction [15,31] or cell shape 
dependent adhesion [29,32], because the basic principles under- 
lying these models are equivalent [12,24]. As a solution to this 
problem, a sufficiently correct complete description of endothelial 
cell behavior should suffice for the emergence of the subsequent 
levels of organization of the system, an approach that requires that 
the system has been experimentally characterized at all levels of 
organization. 

The role of cell traction and ECM mechanics during in vitro 
angiogenesis have been characterized experimentally particularly 
well, making it a good starting point for such a multiscale 
approach. Endothelial cells apply traction forces on the extracel- 
lular matrix, as demonstrated by a variety of techniques, e.g., 
wrinkle formation on elastic substrates [9], force-generation on 
micropillar substrates [33], and traction force microscopy [6,34]. 
Using scanning electron microscopy, Vernon and Sage [9] found 
that ECM ribbons radiate from endothelial cells cultured in 
Matrigel, suggesting that the traction forces locally reorient the 
extracellular matrix. The cellular traction forces produce local 
strains in the matrix, which can affect the motility of nearby cells 
[2] . Thus endothelial cells can both generate, and respond to local 
strains in the extracellular matrix, suggesting a feedback loop that 
may act as a means for mechanical cell-cell communication [2] 
and hence coordinate collective cell behavior. Here, we use a 
hybrid cellular Potts and finite element model to show that a set of 
assumptions mimicking mechanical cell-cell communication via 
the ECM suffices to reproduce observed single cell behavior 
[35,36], pairwise cell interactions [2], and collective cell behavior: 
network formation and sprouting. 

Results 

Response of endothelial cells to static strains in ECM 

First we set out to capture, at a phenomenological level, the 
response of endothelial cells to static strains in the ECM in absence 
of cellular traction forces. When grown on statically, uniaxially 
stretched collagen-enriched scaffolds, murine embryonic heart 
endothelial (H5V) cells orient in the direction of strain, whereas 
cells grown on unstrained scaffolds orient in random directions 
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[37]. Because the collagen fibers make the scaffold stiffen in the 
direction of strain, we hypothesized that the observed alignment of 
cells is due to durotaxis, the propensity of cells to migrate up 
gradients of substrate rigidity [38] and to spread on stiff substrates 
[39,40]. In our model we assumed (a) strain stiffening: a strained 
ECM is stiffer along the strain orientation than perpendicular to it, 
such that (b) due to durotaxis the endothelial cells preferentially 
extend pseudopods along the strain orientation, along which the 
ECM is stiffest, giving cells the most grip. To keep the ECM 
mechanics simulations computationally tractable, we assumed an 
isotropic and linearly elastic ECM. With these assumptions it is not 
possible to model strain stiffening explicitly. We therefore 
mimicked strain stiffening by assuming that stiffness is an 
increasing, linear function of the local strain. 

Durotaxis was modelled as follows, to reflect the observation 
that focal adhesion maturation occurs under the influence of local 
tension [41]: At low local stiffness, we applied standard cellular 
Potts dynamics to mimic the iterative formation and breakdown of 
ECM adhesions, producing "fluctuating" pseudopods. However, if 
the stiffness was enhanced locally, we assumed that the resulting 
tension in the pseudopod led to maturation of the focal adhesion 
[41,42], stabilizing the pseudopod as long as the tension persists. 
To mimic such focal adhesion maturation in the cellular Potts 
model, we increased the probability of extension along the local 
strain orientation, and reduced the probability of retraction (see 
Methods for detail). 

Figure 1 A shows the response of the simulated cells to uniaxial 
stretch along the vertical axis. With increasing values of the 
durotaxis parameter ^durotaxis ( see Eq. 8), the endothelial cells 
elongate more. To test the sensitivity of the durotaxis model for 
lattice effects, we varied the orientation of the applied strain over a 
range [0 — 180] ° and measured the resulting orientation of the 
cells. Figure 1 shows that the average orientation of the cells 
follows the orientation of the stretch isotropically. Thus the 
durotaxis component of our model phenomenologically repro- 
duces published responses of endothelial cells to uniaxial stretch 
[37]. 




' ///III I ' / I \ \ 



Figure 2. Visualization of simulated traction forces [black 
arrows) and resulting matrix strains [blue line segments) 
generated in the proposed hybrid cellular Potts and finite 
element simulation model. 

doi:1 0.1 371/journal.pcbi.l 003774.g002 

Generation of strains in ECM due to cellular traction 

We next attempted to mimic the forces applied by cells onto the 
extracellular matrix, in absence of durotaxis. Traction-force 
microscopy experiments [34,39] show that endothelial cells 
contract and exert tensional forces on the ECM. The forces are 
typically directed inward, towards the center of the cell, and forces 
concentrate at the tips of pseudopods. A recent modeling study by 
Lemmon and Romer [43] found that an accurate prediction of the 
direction and relative magnitudes of these traction forces within 
the cell can be obtained by assuming that each lattice node i 
covered by the cell pulls on every other node the cell covers, j, with 
a force proportional to their distance, djj. Because this model gives 
experimentally plausible predictions for fibroblasts, endothelial 
cells, and keratocytes [43], we adopted it to mimic the cell-shape 
dependent contractile forces that endothelial cells exert onto the 
ECM. Figure 2 shows the contractile forces (black) and resulting 
ECM strains (blue) generated in our model by two adjacent cells. 
The traction forces and ECM strains become largest at the cellular 
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Figure 1. Simulated cellular responses to static strains. Cells do not generate traction forces In this figure. (4) Cell length as a function of the 
durotaxis parameter, laurotaxis/ on a substrate stretched along the vertical axis. (8) Cell orientation as a function of the stretch orientation (simulated 
with Adurotaxis = 10). Error bars show standard deviation for n= 100. Insets show five simulations per value tested. 
doi:1 0.1 371 /journal.pcbi.1 003774.g001 
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"pseudopods", qualitatively agreeing with traction force fields 
reported for endothelial cells [34], 

Mechanical cell-ECM feedback qualitatively reproduces 
effect of substrate stiffness on cell shape and motility 

The two previous sections discussed how the simulated cells can 
respond to and induce strain in the ECM in an experimentally 
plausible way. To test how the simulated cells respond to the 
strains they generate themselves, we studied the behavior of 
simulated, single cells in presence of both the cell traction 
mechanisms and the durotaxis mechanisms. During each time 
step, we used the Lemmon and Romer [43] model to calculate 
traction forces corresponding to current cell positions. Next, we 
started the finite element analysis from an undeformed matrix, 
calculating steady-state strains for the current traction forces. To 
simulate cell movement, which was biased by the local matrix 
strains using the durotaxis mechanism, we then applied one cell 
motility simulation time step, or Monte Carlo step (MCS; the 
MCS is the unit of time of our simulation; see Methods for detail 
and Discussion for an estimate of the real time corresponding to an 
MCS). After running the CPM for one MCS we again relaxed the 
matrix such that the next step started with an undeformed matrix. 
Thus we currently did not consider cell memory of substrate 
strains. 

As Figure 3 and Video SI demonstrate, in this model matrix 
stiffness affects both the morphology and motility of the simulated 
cells. On the most compliant substrate tested (0.5 kPa) the 
simulated cells contract and round up, whereas cells spread 
isotropically on the stiffest substrate tested (32 kPa). Overall, the 
cellular area increases with substrate stiffness (Figure 3 B). On 
matrices of intermediate stiffnesses (around 12 kPa) the cells 
elongate, as reflected by measurements of the cell length (Figure 3 
C) and eccentricity (Figure 3 D) that both have maximum values at 
around 1 2 kPa. Such a biphasic dependence of cellular morphol- 
ogy on the stiffness of the ECM mimics the behavior of endothelial 
cells [39] and cardiac myocytes [36] in matrices of varying 
stiffness. The dependence of cell shapes on substrate stiffnesses is 
due to the transition from fluctuating to adherent pseudopods with 
increasing stiffness. Focal adhesions of cells on soft substrates all 
remain in the "fluctuating" state, irrespective of the local strains. 
On intermediate substrates, some pseudopods, due to increased 
traction, move to an extended state (mimicking a mature focal 
adhesion), generating more traction in this direction. Hence an 
initial stochastic elongation self-enhances in a feedback loop of 
increasing traction and strain stiffening. Such a self-enhancing cell- 
elongation starting from an initial anisotropy in cell spreading has 
previously been suggested by Winer et al [44]. Extensions 
perpendicular to the long axis of an elongated cell do not occur 
since there is insufficient traction and the volume constraint is 
limiting. At matrices of high stiffness all pseudopods attempt to 
extend, mimicking the formation of static focal adhesion, until the 
volume constraint becomes limiting. This makes the cells spread 
more on stiff substrates than on soft substrates, with weaker 
volume constraints (lower values of X ) producing a stronger effect 
of substrate stiffness on cell shape and cell area (Figure SI). 

We also measured the random motility of the cells by 
characterizing their dispersion coefficients, which we derived from 
the mean square displacements of the cells (Figure S2; see section 
Morphometry for detail). The dispersion coefficients show biphasic 
behavior, with the highest motilities occurring at around 12 kPa 
(Figure 3 E). The biphasic dependence of the dispersion to 
substrate stiffness is in accordance with in vitro behavior of 
neutrophils [45] , and smooth muscle cells [46] . Here it is typically 
thought to be due to a balance of adhesion and actin 



polymerization, or due to the interplay between focal adhesion 
dynamics and myosin-based contractility [45]. In our model, the 
effect is more likely due to the appearance of eccentric cell shapes 
at intermediate stiffnesses; as a result, only the tips of the cell 
generate sufficient strain in the matrix to extend pseudopods, 
producing more persistent motion than the round cells at stiff or 
soft substrates. It will be interesting to see if a similar relationship 
between cell shape and cell motility holds in vitro. Thus the model 
rules for cell traction and stretch guidance based on durotaxis and 
strain stiffening suffice to reproduce an experimentally plausible 
cellular response to matrix stiffness. 

Mechanical cell-ECM feedback coordinates behavior of 
adjacent cells 

Strains induced by endothelial cells on a compliant substrate 
with low concentrations of arginine-glycine-aspartic acid(RGD)- 
containing nonapep tides can affect the behavior of adjacent cells 
[2]. On soft substrates (5.5 kPa or below) the cells reduced the 
motility of adjacent cells, whereas on stiff substrates (33 kPa) such 
an effect was not found. On substrates of intermediate stiffness 
(5.5 kPa), adjacent endothelial cells repeatedly attached and 
detached from one another, and cells moved more slowly in close 
vicinity of other cells, than when they were on their own. Because 
the extent to which cells could affect the motility of nearby cells 
depended on matrix compliancy, mechanical traction forces could 
act as a means for cell-cell communication [2]. To test if the simple 
strain-based mechanism represented in our model suffices for 
reproducing such mechanical cell-cell communication, we initiated 
the simulations with pairs of cells placed adjacent to one another at 
a distance of fourteen lattice sites corresponding to a distance of 
35 //m, and ran a series of simulations on substrates of varying 
stiffness (Figure 4 A and Video S2). 

The cells behaved similar to the single cell simulations 
(Figure 3), with little cell-cell interactions at the lower and higher 
stiffness ranges. Consistent with previous observations [2], cell 
pairs on substrates of intermediate stiffness (12 kPa) dispersed 
more slowly than individual cells (paired two-sample £-test at 5000 
MCS, jf?<0.05 for 12 kPa), whereas individual cells and cell pairs 
dispersed at indistinguishable (p<0.05) rates on stiff (14 kPa or 
more) or soft (10 kPa or below) substrates (Figure 4, B-D) and 
Figure S3). 

Also in agreement with the previous, experimental observations 
[2], on a simulated substrate of intermediate stiffness (12 kPa) the cells 
responded to the matrix strains induced by the adjacent cell by 
repeatedly touching each other, and separating again (Figure 4 E). 
The contact duration of cells on soft and stiff substrates, when they get 
close enough to each other, are typically longer than for intermediate 
substrates. This behavior is also similar to observations in vitro [2]. As 
one might expect that strongly adherent cells will not repeatedly 
touch and retract, but rather stay connected upon first contact, we 
investigated the effect of cell adhesion on these parameters (Figure 
S4). Consistent with this intuition, for stronger adhesion, the contact 
count tends to be reduced and the contact durations tend to increase, 
but the overall trend holds: at intermediate matrix stiffnesses we 
continue to observe more frequent cell contacts than for more soft or 
more stiff matrices. Thus the observed pairwise cell behavior is 
primarily driven by durotaxis. 

Mechanical strain can also coordinate the relative orientation of 
cells. Fibroblasts seeded on a compliant gel tend to align in a head- 
to-tail fashion along the orientation of mechanical strain [47]. 
Bischofs and Schwarz [48] proposed a computational model to 
explain this observation. Their model assumes that cells prefer the 
direction of maximal effective stiffness, where the cell has to do the 
least work to build up a force. This work is minimal between two 
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Figure 3. Simulated individual cell responses to mechanical cell-ECM feedback. {A) Single cells on substrates of varying stiffness after 100 
MCS. Line pieces indicate strain magnitude and orientation. (£) cell area {a( a ) ) of cells; (Q cell length (length of major axis if the cell is seen as an 
ellipse) as a function of substrate stiffness (D) cell eccentricity {£, = y/l — b 2 /a 2 , with a and b the lengths of the cell's major and minor semi-axes) 
as a function of stiffness. Mean and standard deviation shown for n= 100 in panels B-D. (£) Dispersion coefficients of individual, simulated cells, 



derived from a linear fit on the mean square displacements (Figure S2); n = 
doi:1 0.1 371 /journal.pcbi.1 003774.g003 



1000. Error bars indicate 95% confidence intervals of linear fits. 



aligned cells, because maximum strain stiffening occurs along the axis 
of contraction. Interestingly, visualization of our model results 
(Figure 1 Q suggested similar head-to-tail alignment of our model 
cells at around 1 2 kPa. To quantify cell alignment in our simulations, 
we measured the angle a between the lines l\ and I2, defining the long 
axes of the cells (Figure 4 F). We classified the angles as acute 
(a < n/ 2; i.e. no alignment) or obtuse (a > n/2; alignment). At matrix 
stiffnesses up to around 10 kPa, about one fourth of the angles a were 
obtuse, corresponding to the expected value for uncorrelated cell 
orientations. However, at 12 kPa and 14 kPa significantly more than 
a fourth of the angles a between the cell axes were obtuse (55/100 for 
12 kPa, p<\ xlO" 8 and 52/100 for 14 kPa, p<\ xlO" 8 , binomial 
test), and for substrate compliancies of 8 to 16 kPa significantly more 
of the angles a were obtuse than for 4 kPa (p<0.0 \ for 8 kPa, andjf?< 
1x10 12 for 10 kPa to 16 kPa; two-tailed Welch's t-test), suggesting 
that the mechanical coupling represented in our model causes cells to 
align in a head-to-tail fashion. 

Mechanical cell-cell communication drives biologically- 
realistic collective cell behavior 

After observing that the local, mechanical cell-ECM interactions 
assumed in our model sufficed for correctly reproducing many aspects 
of the behavior of individual endothelial cells on compliant matrices 
and of the mechanical communication of pairs of endothelial cells on 
compliant matrices, we asked what collective cell behavior the 
mechanical cell-cell coordination produced. When seeded subcon- 
fluently onto a compliant matrix (e.g., Matrigel), endothelial cells tend 
to organize into polygonal, vascular-like networks [5,6,49,50]. To 
mimic such endothelial cell cultures, we initialized our simulations with 
(approximately) 450 cells uniformly distributed over a lattice of 
300x300 pixels (0.75x0. 75mm 2 ), corresponding to a cell density of 
800 endothelial cells per mm 2 . In accordance with experimental 
observations on gels with low concentrations of collagen [6] or RGD- 
peptides [2], after 3000 MCS networks had not formed on soft 
matrices (0.5-4 kPa) or on stiff matrices (16-32 kPa) (Figure 5 A): The 
cells tended to form small clusters (Figure 5 A). Interestingly, on 
matrices of intermediate stiffness after around 300 MCS the cells 



organized into chains (8 kPa) or network-like structures (10 kPa and 
12 kPa) similar to vascular network-like structures observed in 
endothelial cell cultures [5,6,49,50]. The optimal stiffness (~10kPa) 
for network formation is slightly lower than the stiffness of the substrate 
(~12kPa) on which single cells elongate the most (Figure 3 A). In 
comparison with a single cell, the collective pulling of a cell colony 
creates larger strains in the substrate. Consequently, the strain 
threshold inducing cell elongation is crossed at smaller substrate 
stiffness. 

Figure 5 B and Video S3 show a time-lapse of the development 
of a network configuration on a substrate of lOkPa. The cells 
organized into a network structure within a few hundred MCS. 
The network was dynamically stable, with minor remodeling 
events taking place, including closure and bridging of lacunae. 
Figure 5 C shows such a bridging event in detail. In an existing 
lacuna (1800 MCS) stretch lines bridged the lacuna, and 
connected two groups of cells penetrating the lacuna (1980 
MCS). The cells preferentially followed the path formed by these 
stretch lines (2150 MCS) and reached the other side of the lacuna 
by 2400 MCS. Such bridging events visually resemble sprouting in 
bovine endothelial cell cultures on compliant matrices (Figure 5 D, 
Video S4, and [6]). To stay close to the experimental conditions 
used for the observations of pairwise endothelial cell-cell interac- 
tion on compliant substrates [2] that we compared the simulations 
of pairwise interactions with, in this experiment we used a 2.5 kPa 
gel functionalized with 5 jug/ ml RGD peptide - a stiffness at which 
no network-formation is found in our simulations. Although we 
thus do not yet reach full quantitative agreement between model 
and experiment, note that network formation occurs at substrate 
stiffness of lOkPa on polyacrylamide matrices enriched with a low 
(1 jug/ ml) concentration of collagen [6]. 

We next asked if the mechanical model could also reproduce 
sprouting from endothelial spheroids [7,8]. Video S5 and Figure 6 
shows the results of simulations initiated with a two-dimensional 
spheroid of cells after 3000 MCS. On soft (0.5-8 kPa) and on stiff 
(32 kPa) matrices the spheroids stayed intact over the time course 
of the simulation. On matrices of intermediary stiffness (10— 
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Figure 4. Simulated cell-cell interactions on substrates of varying stiffnesses. 04) Visualization of cell shapes and substrate strains in 
absence of external strain. Line pieces indicate strain magnitude and orientation. {B-D) Mean square displacement of individual cells {blue errorbars) 
and cell pairs {red errorbars) on simulated substrates. {B) 4 kPa; (Q 12 kPa; (D) 32 kPa. Error bars in panels B to D indicate standard deviation for 
n = 100. (£) Number of cell-cell contacts made over 500 MCS between two simulated cells initiated at a distance of fourteen lattice sites from each 
other. Error bars show standard deviation over n= 100 simulations (F) Quantification of head-to-tail alignment of cells. An obtuse angle between the 
two cells' long axes indicates that cells are oriented head-to-tail. Plotted is the fraction of Monte Carlo steps over MCS 20-500 that the two cells are 
aligned head-to-tail. Shown are means and standard deviations over 100 independent simulations on a field of 0.25x0.25 mm 2 (100x100 pixels). 
Insets: examples of acute (left) and obtuse (right) cell configurations. 
doi:1 0.1 371 /journal.pcbi.1 003774.g004 



12 kPa) the spheroids formed distinct sprouts, visually resembling 
the formation of sprouts in in vitro endothelial spheroids [7,8]. On 
the 14 kPa and 16 kPa matrices the cells migrated away from the 
spheroid, with some cell alignment still visible for the 14 kPa 
matrices. Observation of a sprout protruding from a spheroid at 
1 0 kPa suggests that a new sprout starts when one of the cells at 
the edge of the cluster protrudes and increases the strain in front of 
it. In a positive feedback loop via an increase in perceived stiffness 
the strain guides the protruding cell forward. The strain in its wake 
then guides the other cells along (Figure 6 C). 

Discussion 

In this paper we introduced a computational model of the in 
vitro collective behavior of endothelial cells seeded on compliant 



substrates. The model is based on the experimentally supported 
assumptions that (a) endothelial cells generate mechanical strains 
in the substrate [34,43], (b) they perceive a stiffening of the 
substate along the strain orientation, and (c) they extend 
preferentially on stiffer substrate [37]. Thus, in short, the 
assumptions are: cell traction, strain stiffening, and duro taxis. 
The model simulations showed that these assumptions suffice to 
reproduce, in silico, experimentally observed behavior of endo- 
thelial cells at three higher level spatial scales: the single cell level, 
cell pairs, and the collective behavior of endothelial cells. In 
accordance with experimental observation [36,39], the simulated 
cells spread out on stiff matrices, they contracted on soft matrices, 
and elongated on matrices of intermediate stiffness (Figure 3). The 
same assumptions also sufficed to reproduce experimentally 
observed pairwise cell-cell coordination. On matrices of interme- 
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Figure 5. Simulated network formation assay. (/A) Simulated collective cell behavior on substrates of varying stiffness, with a uniformly 
distributed initiated configuration of cells. (£) Time lapse showing the development of a polygonal network on a lOkPa substrate (time in MCS). 
Panels A and B represent a 0.75x0.75 mm 2 area (300 x 300 pixels) initiated with 450 cells. (Q Close-up of simulated network formation on a 10 kPa 
substrate, showing the reconnection of two sprouts. Time in MCS. (D) Time lapse imaging of bovine aortic endothelial cells seeded onto a 2.5 kPa 
polyacrylamide gel functionalized with RGD-peptide. Arrows indicate cells that join together and elongate into a network. Time scale is in hours. Scale 
bar is 50 ^m. 

doi:1 0.1 371 /journal.pcbi.1 003774.g005 



diate stiffness, endothelial cells slowed down each other (Figure 4 
B) and repeatedly touched and retracted from each other (Figure 4 
E and Video S2), in agreement with in vitro observations of bovine 
aortic endothelial cells on acrylamide gels [2]. Also, in agreement 
with experimental observations of fibroblasts on compliant 



substrates [47] and previous model studies [48] the cells 
repositioned into an aligned, head-to-tail orientation (Figure 4 
F). The model simulations further suggest that these pairwise cell- 
cell interactions suffice for vascular-like network formation in vitro 
(Figure 5) and sprouting of endothelial spheroids (Figure 6). 
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Figure 6. Simulated spheroid assay. f/\J Collective behavior in a simulation initiated with a two-dimensional "spheroid" of cells, on substrates of 
varying stiffness. (B) Time lapse showing a sprouting spheroid on a lOkPa substrate. Time in MCS. Panels A and B represent a 0.75x0.75 mm 2 area 
(300x300 pixels) initiated with a spheroid consisting of 113 cells; (C) Close-up of sprouting on a 10 kPa substrate. Time in MCS. Black line pieces 
indicate strain magnitude and orientation. 
doi:1 0.1 371 /journal.pcbi.1 003774.g006 



The correlation between pairwise cell-cell interactions and 
collective cell behavior observed in our computational model 
parallels observations in vitro. Cells elongate due to positive 
feedback between stretch-guided extension and cell traction, as 
previously suggested by Winer et al. [44] . Elongated and spindle- 
shaped cells are considered indicative of future cell network 
assembly [6]. Our model suggests that the elongated cell shapes 
produce oriented strains in the matrix, via which cells sense one 
another at a distance. In this way new connections are 
continuously formed over "strain bridges" (see, e.g., Figure 5 
C,D and Video S4), while other cellular connections break, 
producing dynamically stable networks as illustrated in Video S3. 
Such dynamic network restructuring was also observed during 
early embryonic development of the quail embryo [51] and in 
bovine aortic endothelial cell cultures (Figure 5 D and [6]), but not 



in human umbilical vein endothelial cell cultures [23,50]. Also in 
agreement with experimental results, the collective behavior 
predicted by our model strongly depends on substrate stiffness. 
The strongest interaction between cell pairs is found on substrates 
of intermediate stiffness, enabling network formation [2], whereas 
network assembly does not occur on stiffer or on softer 
substrates [6] . 

These agreements with experimental results are encouraging, 
but our model also lacks a number of properties of in vitro 
angiogenesis that pinpoint key components still missing from our 
description. We compared the simulation of pairwise cell-cell 
interactions with previous experiments conducted on polyacryl- 
amide gels, functionalized with RGD ligands [2], which have 
linear elastic behavior for small deformations [52-54]. Strain- 
stiffening of polyacrylamide gels has been reported for deforma- 
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tions over 2 jUm [55]. Thus with pixels in our model measuring 
2.5 //mx2.5 //m, strain-stiffening seems a reasonable assumption. 
Nevertheless, a possible alternative interpretation of the cell pair 
simulations is that the increased tension generated in pseudopods 
pulling on the matrix leads to a higher probability of focal 
adhesion maturation [4 1 ,42] . A further issue is that in our 
simulations, single cells dispersed somewhat more quickly on soft 
gels than on stiff gels (Figure 3 E and Figure S2). This model 
behavior contradicts experimental observations that endothelial 
cells move fastest on stiff substrates [2]. Another open issue 
concerns the time scales of our simulations. In the present paper 
time we use the Monte Carlo step as a (computational) unit of 
time. To estimate the actual time corresponding to 1 MCS, we 
scale the single cell dispersion coefficients shown in Figure 3 E to 
experimental dispersion coefficients of bovine endothelial cells on 
compliant substrates in vitro [2]. Reported dispersion coefficients 
of endothelial cells range from around 1 /im 2 /min (on substrates 
of 500 Pa) to around 10 /rni 2 /min (on substrates of 5500 Pa) (as 
derived from the MSDs in Figure 3a,c in [2] and based on 
MSD(7) = 4Z)/; cf. Eq. 13). The dispersion coefficients of single 
cells in our simulations are in the range of 0.03 — 0.08 /im 2 /MCS 
(Figure 3), assuming pixels of 2.5 x 2.5 /mi 2 . Thus, based on fitting 
of single cell dispersion rates, the estimated length of 1 MCS is 0.5 
to 3 seconds. The typical time scale of a vascular network 
formation simulation is around 3000 MCS (Figure 5), i.e., 
12.5 min to 2.5 hr for these time scale estimates. In experiments, 
network formation takes longer, around 24 hr. Thus in our 
current model the time scales of cell dispersion and network 
formation do not match exactly. A possible reason of this 
discrepancy is the short persistent length of cell motility in 
standard cellular Potts models. To better match the time scales of 
single cells and collective cell behavior in our model, in our future 
work we will increase the persistence length of the endothelial cells 
by using the available cellular Potts methodology [56-58], or 
model the subcellular mechanisms of cell motility in more detail, 
e.g. by including mean-field models of actin polymerization 
[59,60]. A further open issue is the interaction between substrate 
mechanics and cell-substrate adhesivity. Although the model 
correctly predicts the absence of network formation on stiff 
substrates, it cannot yet explain the observation that reducing the 
substrate adhesivity of the endothelial cells rescues network 
formation on stiff substrates [6]. On compliant gels endothelial 
cells must secrete fibronectin to form stable networks, whereas 
fibronectin polymerization inhibitors elicit spindle-like cellular 
phenotypes associated with network formation on stiff matrices, 
under conditions where networks do not normally form [6]. To 
explain these observations, straightforward future extensions of the 
model will include a more detailed description of cell-substrate 
adhesion, combined with models of ECM secretion and proteolysis 
[13,25,27,61]. 

The current model also assumes a uniform density (i.e., the 
infinitesimal strain assumption) and thickness of the extracellular 
matrix, whereas under some culture conditions the endothelial 
cells have been reported to pull the extracellular matrix 
underneath them [62], producing gradient in matrix density 
and/ or thickness. To describe the role of viscous deformations of 
the extracellular matrix in morphogenesis, Oster and Murray 
[17,18] developed a continuum mechanical model of pattern 
formation in mesenchymal tissues. Their model assumed (a) that 
cells exert contractile forces onto the surrounding extracellular 
matrix, that will (b) locally deform the ECM, resulting in passive 
displacements of cells along with the ECM, and (c) produce density 
gradients in the ECM along which cells move actively due to 
haptotaxis. These mechanisms together produce periodic cell 



density patterns. Manoussaki et al. [15] and Namy et al.[19] 
applied this work to investigate mechanical cell-ECM interactions 
during angiogenesis, and demonstrated that the mechanism can 
produce vascular-like network patterns. In their model they also 
included an anisotropic diffusion term to simulate preferential 
movement along the local strain-direction, but the term was 
neither necessary nor sufficient for network formation. This 
finding contradicts our model in which strain-induced sprouting is 
the driving force of network formation and sprouting. The two 
models represent the two extremes of network formation on visco- 
elastic matrices. Here, the Manoussaki et al. [15] and Namy et al. 
[19] models represent patterning on viscous matrices, in which 
cellular traction forces pull the matrix together while inducing little 
strain or stress. Our model would represent elastic materials, in 
which pulling forces induce local strains. Future extensions of the 
model will include matrix remodelling (e.g., by assuming a matrix 
thickness field) allowing us to study the full range of viscoelastic 
matrices. 

Apart from these biological issues, we made several mathemat- 
ical simplifications that we will improve upon in future models of 
cell-ECM interactions. In the current model, for mathematical 
simplicity, we assumed that after each Monte Carlo step the 
matrix was undeformed again. Thus we currently did not consider 
cell memory of substrate strains. Further developments of the 
model presented here will improve on this issue, because actin 
filament dynamics are typically influenced by the past evolution of 
substrate deformations, e.g., due to reorientation of matrix fibers 
[62]. For computational efficiency, we assumed linearly elastic 
materials and infinitesimal strain in the finite element simulations, 
and mimicked durotaxis via a perceived strain-stiffening (Eq. 9) 
where cells perceive increased ECM stiffness due to local strain. In 
our ongoing work we are interfacing the open source package 
FEBio (http://febio.org) with the cellular Potts package Compu- 
Cell3D (http://compucell3D.org). This will allow us to run our 
model with any ECM material available to users of FEBio, 
including strain-stiffening materials. Using an actual strain 
stiffening material may lead to longer-range interactions between 
cells, because locally stiffer regions may channel the stress between 
the cells [63]. A further technical limitation of our model is that we 
currently only run two-dimensional simulations, representing cells 
moving on top of a two-dimensional culture system. The ongoing 
interfacing of FEBio and CompuCell3D will pave the way for 
modeling cell-ECM interactions in three-dimensional tissue 
cultures. We also plan to model fibrous extracellular matrix 
materials in more detail. 

A quite puzzling aspect of vascular network formation and 
spheroid sprouting is that so many alternative, often equally 
plausible computational models can explain it (reviewed in [12]). 
Including the present model, there are at least three alternative 
computational models based on mechanical cell-ECM interactions 
[15,16,19,31,64], a series of models assuming chemoattraction 
between endothelial cells [20,21,23,24,65,66] and extensions 
thereof [25,27,67], and models explaining network formation in 
absence of chemical or mechanical fields [29,30,32]. Each of the 
models explains one aspect of vascular network formation or a 
response to an experimental treatment that the other models 
cannot explain, e.g. the relation between spindle-shaped cell 
phenotypes and network formation [23,30], the requirement of 
VE-cadherin signaling for network formation and sprouting 
[24,29], the binding and release of growth factors from the 
ECM [25,26], the role of mechanical ECM restructuring and 
haptotaxis [15,19,31], the response of vascular networks to toxins 
[27], or the role of intracellular Ca 2+ signaling [57]. Among these 
alternative models, we must now experimentally falsify incorrect 
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mechanisms, and fine-tune and possibly combine the remaining 
models to arrive at a more complete understanding of the 
mechanisms of angiogenesis. To this end, we are currently 
quantitatively comparing the kinetics of patterns produced by 
chemotaxis-based, traction-based, and cell-elongation based mod- 
els with the kinetics of in vitro networks [23,50]. The resulting, 
more complete model would likely contain aspects of each of the 
available computational models and assist in explaining the 
conflicting results obtained from the available experimental 
systems, culture conditions, and in silico models of angiogenesis. 

Methods 

To model the biomechanical interactions between endothelial 
cells and compliant matrices, we developed a hybrid of the cellular 
Potts model (CPM) [68,69] to represent the stochastic motility of 
the endothelial cells, and a mechanical model based on the finite 
element method (FEM) [70] of the compliant extracellular matrix. 
Related CPM-FEM models were proposed for the simulation of 
load-induced bone remodeling [71,72], and recently a related 
approach was proposed in a model study of cell alignment [73]. A 
documented simulation code is provided as part of the Supporting 
Information (Supporting Text SI and Code SI) and a detailed list 
of parameter values is given in Table S 1 . 

Cellular Potts model 

The CPM represents cells on a regular square lattice, with one 
biological cell covering a cluster of connected lattice sites. To 
mimic random cell motility, the CPM iteratively expands and 
retracts the boundaries of the cells, depending on the passive forces 
acting on them and on the active forces exerted by the cells 
themselves. These are summarized in a balance of forces, 
represented by the Hamiltonian, 



of the cells to their environment, the copy step is always accepted if 
it decreases the Hamiltonian. To account for the active random 
motility of biological cells, we allow for energetically unfavorable 
cell moves, by accepting copies that increase the Hamiltonian with 
Boltzmann probability, 



P(AH) = 



1 

e -AH/T 



if AH<0 
if AH>0, 



(2) 



where AH is the change in H if the copying were to occur, and 
T>0 parameterizes the intrinsic cell motility. It represents the 
extent to which the active cell motility can overcome the reactive 
forces (e.g. volume constraint or adhesions) in the environment. 
We assume that all cells keep the same motility and thus set T to 
be constant throughout the simulations. During one Monte Carlo 
step (MCS), we perform n copy attempts, with n equal to the 
number of sites in the lattice. To prevent cells from splitting up 
into two or more disconnected patches, we use a connectivity 
constraint that rejects a spin flip (j(x')^-x if it would break apart 
the retracting cell o(x). 

Model of compliant substrate based on finite element 
method 

A two-dimensional model describes the compliant substrate on 
which the cells move. Deformations are calculated using the finite 
element method (FEM; reviewed in [70]). The FEM represents the 
substrate as a lattice of finite elements, e, with each element 
corresponding to a pixel of the CPM. To obtain the finite element 
equations, the weak formulation (associated with the total potential 
energy) of the governing equations of the displacement u of the 
substrate is set up, in order to obtain the finite element equations, 



H = 

aecells V / (x,x f ) 

The first term is an (approximate) volume constraint, with a(o) 
the actual volume of the cells, A(a), a resting volume, and X an 
elasticity parameter that regulates the permitted fluctuation 
around the resting volume. In contrast with the original 
formulation of the CPM [68], the deviation of the cell from its 
target volume is taken relative to the target volume, by analogy 
with the (non-dimensional) engineering strain. Alternative, similar 
volume constraints can be chosen [67]. We use a value A(<j) = 50 
for all cells; the medium does not have a volume constraint. The 
second term represents cell-cell and cell-medium adhesion, where 
J(g(x) 9 g(x?))>0 is the contact cost between two neighboring 
pixels, and 3, the Kronecker delta. Throughout the manuscript we 
use neutral cell-cell adhesion settings; J((t(x),<j(x')) = 2.5 at cell- 
cell interfaces, and J((t(x),0) = 1 .25 at cell-medium interfaces, 
with <j(x)>0 and cr(x / )>0. In other words, cells have no 
preference for adhering to other cells or the medium. For these 
neutral cell adhesion parameter settings, cells will still adhere 
weakly to one another (a remedy for this effect was proposed in 
[74]). Additional terms in the Hamiltonian represent the cells' 
responses to ECM mechanics, and will be described in more detail 
below. 

The CPM iteratively selects a random lattice site x' and 
attempts to copy its state, o(x'), into a randomly selected adjacent 
lattice site x. To reflect the physical, "passive" behavioral response 



Ku=f, (3) 

with stiffness matrix K, displacement u, and forces f. The vector 

u= [u Xl ,u yi ,u X2 ,...u Xn ,Uy n ] T contains the displacements of all nodes, 
which are the unknowns that the FEM calculates based on the 
active forces exerted onto the material, presented in /. In this 
paper / only consists of traction forces that the cells apply onto the 
ECM, unless stated otherwise. In a two-dimensional analysis the 
forces / are divided by the thickness they are working on. For this 
we assume an effective substrate thickness t= 10 jim. We impose 
boundary conditions of w = 0 at the boundary of the CPM grid, 
this means that the substrate is fixed along the boundaries. 

To a first approximation, in this work we consider an isotropic, 
uniform, linearly elastic substrate [48,75] and we apply infinites- 
imal strain theory: We assume that material properties, including 
local density and stiffness are unchanged by deformations. The 
global stiffness matrix K is assembled from the element stiffness 
matrices K (see Supporting Text SI and [70]), which describe the 
relation between nodes of each element, 

K = [ B T DBdQ e . (4) 

where B — the conventional strain-displacement matrix for a four- 
noded quadrilateral element (see Supporting Text SI and [70]) — 
relates the node displacements u e to the local strains, as, 



PLOS Computational Biology | www.ploscompbiol.org 



10 



August 2014 | Volume 10 | Issue 8 | e1003774 



Modeling Mechanical Cell-Matrix Feedback 



e = Bu 0 . 



(5) 



The strain vector £ is a column notation of the strain tensor e 
and D is the material property matrix. Assuming plane stress 
conditions, 



is assumed to be undeformed again, for the sake of simplicity. 
Thus, the stiffness matrix, K, is constant in time. 

We assume durotaxis, i.e., the CPM cells preferentially extend 
pseudopods on matrices of higher stiffness (e.g., because of strain 
stiffening). By analogy with chemotaxis algorithms [79] at the time 
of copying we add the following durotaxis term to AH in response 
to the strain- and orientation-dependent ECM stiffness E, 



D- 



1 v 
V 1 



(6) 



0 0 i(l-v), y 



where E is the material's Young's modulus, and v is Poisson's 
ratio. Throughout this study, we use a Poisson's ratio v = 0.45 and 
Young's moduli ranging from is =0.5 kPa to E = 32 kPa, which 
are plausible values for most cell culture substrates [48,53,76]. For 
more details of the derivation of Eq. 3, and the entries in B, see 
Supporting Text SI and [70]. 

As a reference configuration for the displacements we used an 
unstretched substrate, w = 0. Thus, after each Monte Carlo step 
(during which the cells positions and shapes have changed) the 
substrate is assumed to be undeformed, such that the stiffness 
matrix, K, is constant in time. This prevents expensive calculations 
that would be necessary for recalculating K in each iteration. 
Although the previous displacements do not influence the new 
deformation of the substrate, they are used as an initial guess for 
solving Ku = f, in order to reduce the number of iterations 
necessary to converge to the FEM solution. 

Mechanical cell-substrate coupling 

To simulate cell-substrate feedback we alternate the cellular 
Potts model (CPM) steps with a simulation of the substrate 
deformations using the finite element method. We assume that 
cells apply a cell-shape dependent traction on the ECM and the 
cells respond to the resulting ECM strains by adjusting their cell 
shape. Using the CPM grid as the finite element mesh, the pixels 
of the CPM become four-node square elements in the FE-mesh. 
Adopting the model by Lemmon & Romer [43], we assume that 
each node i covered by a CPM cell pulls on all other nodes j in the 
same cell, at a force proportional to distance dij. The resultant 
force E; on node i then becomes, 



Ei = fi^2dij, 



(7) 



where Ax is the lattice spacing and fi gives the tension per unit 
length. This parameter has been scaled to /i = 0.01 nN//mi, such 
that the total cell traction corresponds to experimentally reported 
values [77]. The resultant forces point towards the cell centroid, 
and are proportional to the distance from it (Figure 2). In this way 
a CPM configuration yields a traction force F_, which are collected 
in the forces / for the finite element calculation. To calculate the 
resulting ECM strains, we solve Ku=f for the node displacements 
u with a preconditioned conjugate gradient (PCG) solver [78], and 
derive the local strains using Eq. 5. The reference configuration for 
the displacements is an unstretched substrate, w = 0. After a 
sufficiently accurate solution for the FEM equations has been 
obtained by the PCG solver, we run a Monte Carlo step of the 
CPM. After each MCS, which changes cell positions, the substrate 



-g(x,x')^durotaxis (h(E(e\))(y\ 'V m f + h(E(£ 2 ))(v 2 ' v m ) 2 ) , 



(8) 



with g(x,x') = 1 for extensions and g(x,x') = — 1 for retractions, 

^durotaxis is a parameter, v m =x — x f , a unit vector giving the copy 
direction, and e\ and £2, and Vi and v 2 eigenvalues and eigenvectors 
of £ representing the principal strains and strain orientation. We use 
the strain e(x) in the target pixel when considering an extension, 
and for retractions we use the strain in the source pixel, e(x r ). Thus 
we consider the strain in the ECM adjacent to the pseudopod. The 
sigmoid h{E) = 1/(1 + exp ( — fi(E — Eo))), with threshold stiffness 
Eq, and the steepness of the sigmoid, mimics maturation of focal 
adhesions under the influence of tension [41]. The tension in focal 
adhesions will increase with higher local matrix stiffness, E, because 
the matrix will deform less easily. The sigmoid function starts at 
zero, goes up when there is sufficient stiffness, and eventually 
reaches a maximum. This means that a certain level of stiffness is 
needed to cause a cell to spread. Alternative forms of h(E) can be 
used: For an overview see Figure S5. Due to limitations of our 
current finite element code and for reasons of computational 
efficiency, we assumed a linearly elastic, isotropic material in the 
FEM, thus precluding explicit strain stiffening effects in the FEM 
calculations. Instead, we implemented the effect of strain-stiffening 
in the cell response, where cells perceive increased ECM stiffness as 
a function of the principal strains £\ and £2, 



E(e) = E 0 (l+(e/e st )h> 0 ) 



(9) 



where Eq sets a base stiffness for the substrate, and £ st is a 
stiffening parameter. The indicator function l e >o = {1,£>0; 0, 
£<0} indicates that strain stiffening of the substrate only occurs 
for substrate extensions (e>0); compression (e<0) does not stiffen 
or soften the substrate. 

Morphometry 

To characterize the random motility of single cells and cell 
pairs, we measured the cells' mean square displacement, 



MSD(0 = <(C(5,0 - C(Sfi)f), 



(10) 



with C(S,t), the centroid of cell S at Monte Carlo step ("time") t, 
given by 



C(S,t) = 



1 



(ii) 



with C(S,t), the set of coordinates of the lattice sites comprising 
cell S at MCS L 
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C(S,t) = {x : xeZ 2 A a(x,t) = S}, (12) 

and x= {x\,X2}. The MSD is a reliable measure of random 
motility [80] and it can be directly compared with experimental 
data (e.g., [2]). 

The dispersion coefficient, defined as 

D= limi<(C(S,0-C(5,0)) 2 >, (13) 

t->CTj 4t 



is derived from the slope of the MSD, and is used as a measure of 
the motility of random walkers. The length, orientation and 
eccentricity of cells were estimated from the inertia tensors I(S) of 
the cells, defined as [81], 




Efccra C2(5))2 _ " Efccra (XI " C ^ (X2 - C2(5)) \ (1 4) 



Assuming cells are approximately ellipse-shaped, the length of 
cell a is approximated as 1(a) = 4^/e 2 (I(S)) / \ C(S) | , with e 2 (I(a)) 
the largest eigenvalue of I(S). The eccentricity of a cell is defined 
using the eigenvalues of the inertia tensor 1(a) as 

« ff ) = ^l-(^||) 2 , where e ,(/(S))< e2 (/(S)) are the 

eigenvalues of I(S). An eccentricity close to zero corresponds to 
roughly circular cells and cells with an eccentricity close to unity are 
more elongated. The orientation of the cell is given by the 
eigenvectors of the inertia tensor I(S). 

Endothelial cell culture 

Bovine aortic endothelial cells (BAECs) (VEC Technologies, 
Rensselaer, NY) were cultured through passage 12. Cells were 
kept at 37°C and 5% CO2 and fed every other day with 
Medium 199 (Invitrogen, Carlsbad, CA) supplemented with 
10% Fetal Clone III (HyClone, Logan, UT), 1% MEM amino 
acids (Invitrogen), 1% MEM vitamins (Medtech, Manassas, 
VA), and 1% penicillin-streptomyocin (Invitrogen). Polyacryl- 
amide hydrogels were synthesized as previously described [6]. 
Briefly, a gel mixture was prepared from MilliQ, water, 
HEPES, TEMED (Bio-Rad, Hercules, CA) and a 5%:0.1% 
ratio of acrylamide to bis-acrylamide (Bio-Rad) to generate 
substrates with a Young's modulus of 2,500 Pascals. Polymer- 
ization was initiated by the addition of N-6-((acryloyl)amido)- 
hexanoic acid (synthesized according to Pless et al. [82]) and 
ammonium persulfate (Bio-Rad) into the gel mixture. Follow- 
ing polymerization, gels were incubated with 5 fig/ml RGD 
peptide (GCGYGRGDSPG) (Genscript), followed by ethanol- 
amine (Sigma). Gels were stored in PBS overnight. Hydrogels 
were sterilized with ultraviolet light before cell culture. A T-75 
flask with a confluent BAEC monolayer was seeded onto the 
hydrogels at 350,000 cells per gel (approximately 1,375 cells 
per mm 2 ). The gels were maintained at 37°C and 5% CO2 for 
three days prior to imaging. After replenishing with fresh 
complete media, the cells on hydrogels were visualized with a 
Zeiss Axio Observer. Zl inverted spinning disc microscope with 
a Hamamatsu ORCA-R 2 digital camera. Images were 
captured every 30 minutes for 24 hours. 



Supporting Information 

Figure SI Simulated responses of individual cells to 
mechanical cell-ECM feedback as a function of the 
values of the volume restriction, 1. Columns: area (left), cell 
length (middle) and eccentricity (right). Mean and standard 
deviation shown for 72=100 after 500 MCS on simulated 
substrates of stiffness varying from 0.5 kPa to 32 kPa. 
(PDF) 

Figure S2 Mean square displacements of individual 
cells on simulated substrates of stiffness varying from 
0.5 kPa to 32 kPa. Mean square displacement shown over 
w=1000 cells. 
(PDF) 

Figure S3 Mean square displacement of individual cells 
(blue errorbars) and cell pairs (red errorbars) on 
simulated substrates of stiffness varying from 0.5 kPa 
to 32 kPa. Error bars indicate standard deviation for n = 100. 
(PDF) 

Figure S4 Number of cell-cell contacts made over 500 
MCS (left column) and contact duration (right column) 
over 500 MCS between two simulated cells initiated at a 
distance of fourteen lattice sites from each other on 
simulated substrates of stiffness varying from 0.5 kPa to 
32 kPa, for intercellular contact energies varying from 
J(a(x),(T(x')) = 0.5 (adhesive cells) to J(o(x),g(x')) = 4 (re- 
pulsive cells), with a(x)>\ and a(x')>\\ J(g(x),0) = 1.25 
for all simulations. 
(PDF) 

Figure S5 Effect of form of model function h(E) on cell 
shapes on substrates of different stiffnesses. (A) Standard, 
sigmoid function, as used in main text, h(E) = 
1/(1+ exp( £ 0 ))) with £ 0 = 15000, p = 0.0005, and 
^durotaxis = 10. {B) Saturated function, h(E) = (E/Eq)/(1 + E/E 0 ), 
with ^0 = 15000 and Adurotaxis = 25. (C) Piecewise linear function, 
h(E) = {E/oc, E<E max , E>E mSLX }, with E max = 30000, a = 30000, 
and ^durotaxis = 20. (D) Gaussian function, h(E) = 
Qxp(-(E-E 0 ) 2 /(2y 2 )), with £ 0 = 15000 and y = 2000, 
Idurotaxis = 10. Insets show typical cell shape for regions indicated 
by red bars. 
(PDF) 

Protocol SI C and Matlab source code used for the 
simulations. 

(ZIP) 

Table SI Parameter settings of the simulation model. 

(PDF) 

Text SI Documentation of C and Matlab code used for 
the simulations, including a detailed description of the 
finite-element model. 

(PDF) 

Video SI Behavior in silico of a single cell on substrates 
of 4 kPa, 12 kPa, and 32 kPa, for a duration of 500 MCS 
per simulation. Parameter settings as in Figure 3. 
(MOV) 

Video S2 Pairwise cell-cell interactions in silico on 
substrates of 4 kPa, 12 kPa, and 32 kPa, for a duration 

of 500 MCS per simulation. Parameter settings as in Figure 4. 
(MOV) 

Video S3 Network formation in silico on a substrate of 
lOkPa, for a duration of 3000 MCS. Video represents a 
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0.75 x 0.75 mm 2 area (300 x 300 pixels) initiated with 450 cells. 

Parameter settings are as in Figure 5. 

(MOV) 

Video S4 Network formation of bovine aortic endothe- 
lial cells on a 2.5 kPa polyacrylamide gel functionalized 
with RGD-peptide. Time lapse images were captured in 
30 minute intervals over an 8 hour time period. Image size as in 
Figure 5 D. 
(MOV) 

Video S5 Sprouting in silico from a spheroid on a 
substrate of lOkPa, for a duration of 3000 MCS. Video 
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